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Abstract: Over the past few decades, general relativity and the concordance ACDM model 
have been successfully tested using several different astrophysical and cosmological probes 
based on large datasets (precision cosmology). Despite their successes, some shortcomings 
emerge due to the fact that general relativity should be revised at infrared and ultraviolet 
limits and to the fact that the fundamental nature of dark matter and dark energy is still a 
puzzle to be solved. In this perspective, f{R) gravity has been extensively investigated, 
being the most straightforward way to modify general relativity and to overcame some 
of the above shortcomings. In this paper, we review various aspects of f{R) gravity at 
extragalactic and cosmological levels. In particular, we consider a cluster of galaxies, 
cosmological perturbations and N-body simulations, focusing on those models that satisfy 
both cosmological and local gravity constraints. The perspective is that some classes of f{R) 
models can be consistently constrained by the large-scale structure. 

Keywords: f{R) gravity; cluster of galaxies; N-body simulations; cosmological 

perturbations; growth factor; Cosmic Microwave Background 
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1. Introduction 

The concordance ACDM cosmological model is based on Einstein’s general relativity (GR), 
the standard model of particles with the inclusion of two new ingredients, which are the cosmological 
constant A and dark matter. It provides a theoretical framework to describe the formation and the 
evolution of the structures in the Universe. It has been successfully tested using many different datasets, 
such as the Supemovae Type la (SNela) [1-4], the matter power spectrum from the two-degree field 
(2dF) survey of galaxies [5,6] and from the Sloan Digital Sky Survey (SDSS) [7], the temperature 
fluctuations of the Cosmic Microwave Background (CMB) radiation [8-18], the baryonic acoustic 
oscillations (BAG) [19], and so on. Despite its successes, the following questions are still need to be 
answered: Is GR sufficient to explain all of the gravitational phenomena from collapsed objects to the 
evolution of the Universe and the formation of the structures? Does the theory need to be changed, or at 
least modified? There exist two different approaches to solve the debate: preserve the successes of GR by 
incorporating new particles and/or scalar fields not yet observed and improve the accuracy of the data to 
further verify the model; or modify the theory of gravity to make it compatible with quantum mechanics 
and cosmological observations without introducing additional particles and fields. Independently of the 
preferred approach, having alternatives demands testing the model more deeply. In the concordance 
ACDM model, the most important energy density component is the cosmological constant A (or more 
in general, the dark energy (DE)), required to explain the accelerated expansion of the Universe, with 
Ua = 0.691 ± 0.006 [13] in units of the critical density. The second in importance is dark matter (DM), 
with VIdm = 0.259 ± 0.005 [13], needed to explain the emergence of the large-scale structure (ESS). 
However, the lack of a full comprehension of the fundamental nature of those two components, whether 
particles or scalar fields, is completely unknown and has demanded to further test the GR to understand 
if it is the effective theory of gravity. Additionally, the need for requiring two unknown components to 
fully explain the astrophysical and cosmological observations has been interpreted as a breakdown of the 
theory at those scales [20]. 

As an alternative to introducing extra matter components in the stress-energy tensor, it is possible 
to change the geometrical description of the gravitational interaction, adding in the Hilbert-Einstein 
Eagrangian, higher-order curvature invariants, such as R^, R ^R or R O’^R), and 

minimally or non-minimally coupled terms between scalar fields and geometry (such as (j)^R) [20-28]. 
Extended theories of gravity (ETGs) can be classified as: scalar-tensor theories, when the geometry 
is non- minimally coupled to some scalar field; and higher order theories, when the action contains 
derivatives of the metric components of order higher than the second. ETGs differ from GR in many 
concrete aspects and, therefore, are testable with current and forthcoming experiments. The most 
straightforward way to extend the theory of gravitation is replacing the Einstein-Hilbert Eagrangian, 
which is linear in the Ricci scalar, R, with a more general function of the curvature f{R). Since 
the exact functional form of the /(i?)-Eagrangian is unknown, theoretical predictions should be matched 
to all available data (at all scales) to accommodate the observations. f{R) models have been used to 
test: stellar formation and evolution [29-34]; emission of gravitational waves and constraints on the 
massive graviton modes [35-44]; the flat rotation curves of spiral galaxies [45]; the velocity dispersion 
of elliptical galaxies [46]; to describe a cluster of galaxies [47-53]; the structure formation and the 
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evolution of the Universe [14,54-65]. Although all tests provide elear indieations that the /(i?)-gravity 
models eould overeame the shorteomings of the ACDM model, a final and self-eonsistent framework 
has not been reached yet. 

In addition to the simple f{R) model, it is possible to construct much more complex Lagrangians. As 
an example, teleparallel gravity considers Lagrangians that are a function of torsion scalar /(T) [66]. In 
such models, the Ricci scalar is assumed to vanish, and its role is played by torsion. Furthermore, the 
particularity of these /(T) models is that they include torsion as the source of DE without considering 
a cosmological constant, and they give rise to second order field equations that are easier to solve than 
the fourth order ones from f{R) gravity [66-77]. These models have been widely investigated from 
fundamental to cosmological scales [76,78-84]. Recently, many efforts been devoted to investigating 
the f{R,T) models, where the gravitational Lagrangian consists of an arbitrary function of the Ricci 
scalar R and the torsion scalar T [85-88]. Another class of interesting models that have been introduced 
are, for example, f{G), where Q is Gauss-Bonnet topological invariant or combinations of these with the 
Ricci scalar as the f{R, G) [89-96]. Using these models, it is possible to construct viable cosmological 
models, preventing ghost contributions and contributing to the regularization of the gravitational action. 

In this review, we will focus only on f{R) models. We summarize the main achievements, pointing 
out the difficulties and the future perspectives in probing the /(i?)-gravity using LSS datasets. In 
Section 2, we report the general formalism and the main features of /(i?)-gravity models, and also, 
we focus on two models/approaches that have led to several results in the last few years. In Section 3, we 
describe two approaches to probe /(i7)-gravity using a cluster of galaxies. In Section 4, we review the 
most important results obtained using N-body simulations to quantitatively describe the physical effect 
of /(i7)-gravity models on the structure formation. In Section 5, we summarize the main approaches 
to test alternative models using the expansion history of the Universe. In Section 6, we consider the 
constraints obtained using the CMB power spectrum. Finally in Section 7, we give the main conclusions 
and the future perspectives for this field. 

2.f(R) Gravity 

In metric formalism, the simplest four-dimensional action in f{R) gravity is: 

4 R 

-4 = J [f{R) + ^matter] (1) 

where g is the determinant of the metric and ^matter is the standard perfect fluid matter Lagrangian. 
Varying the action Equation (1) with respect to the metric tensor, we obtain the field equations and 
its trace: 

f,RR^iu - ^f{R)9t,u - f^R = ( 2 ) 

IrR - 2fiR) + SDIr = (3) 

Here, f^R = df /dR, □ = g^^V^ is the d’Alembert operator, and ^ 

energy momentum tensor of the matter. We can recast the Equation (2) in Einstein form as follows: 
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G,u = ^ I ij,. !/(«)--«/, 


r] + — Q/iuOf^R > + 
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where we can reinterpret the first term on the right side of the equation as an extra stress-energy tensor 
contribution f{R) gravity is sufficiently general to include all of the basic features of ETGs, but 

at the same time, it can be easily connected to the observations. Let us consider now two classes of f{R) 
models, which are designed to satisfy cosmological and Solar System constraints. 


2.1. Chameleon Models 

One way to modify gravity is to introduce a scalar field (0), which is coupled to the matter components 
of the Universe. The range and the length of interactions mediated by the scalar field depend on 
the density of the environment. In high-density configurations, the effective mass of the scalar field 
is large enough to suppress it and to accommodate all existing constraints from Solar System tests 
of gravity. Since the mass depends on the local matter density, the cosmological evolution of the 
scalar field could be an explanation for the acceleration of the Universe. Many environment-dependent 
screenings have been proposed, such as the chameleon [97], the dilaton [98], the symmetron [99] or the 
Vainshtein [100,101] mechanisms. One of the most tested mechanism is the chameleon; it must satisfy 
the following equation [97]: 

.2. .. , P 


= V.+ 


M, 


-P 


(5) 


PI 


where V is the potential of the scalar field, (3 is the coupling to the matter, p is the matter density and 
Mpi = s/SttG is the Planck mass. The chameleon field leads to the fifth force in the form of: 


m = 


/S 


Ml 


■V0 


PI 


f{R) models show a chameleon mechanism having a fixed value of the coupling (3 = 
have an additional degree of freedom that mediates the fifth force [27]: 


( 6 ) 


. These models 


V 3 Mpi 


(7) 


where 0oo describes the efficiency of the screening mechanism when r —)■ cxo. Thus, the f{R) Lagrangian 
cannot be considered as a completely free function; it must satisfy many different constraints coming 
from both theoretical and phenomenological side. For example, in a region with high curvature, the 
first derivative of the Lagrangian must be f^R > 0, as well as the second derivative (we define: f^RR = 
d^f{R)/dR^) /, RR > 0. These conditions will avoid tachyons and ghost solutions. Furthermore, to 
satisfy. Solar System constraints must be |/,ijo| << 1 in the present Universe. Two well-studied choices 
matching these constraints are the Starobinsky [102] and Hu-Sawicki [103] models. In this review, we 
will focus only on the latter, which is expressed by the following Lagrangian: 

Cl (5f 


f{R) = -m 




( 8 ) 
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from which one can obtain: 

( 2 \ n+1 

m \ 

-r) 

where the Ricci scalar for the flat ACDM model is: 


R 3m^ 



(9) 


( 10 ) 


dyV / o 

Here, the term — matches the standard when ci/cn —)■ 0. Hu and Sawicki [103] have shown 
o O 

that to recover GR results within the Solar System, one must be 1/,^ I < 10“®. 


2.2. Analytical f(R) Gravity Models and Yukawa-Like Gravitational Potentials. 

In this subsection, we will illustrate how an analytic f{R) model expandable in Taylor series 
gives rise to a Yukawa correction to the Newtonian gravitational potential. Phenomenologically, 
this kind of correction behaves as the so-called chameleon mechanism, becoming negligible at small 
scales [104,105]. In the case of the post-Newtonian corrections, such a solution leads to redefining the 
coupling constants in order to fulfill the experimental observations. In order to derive the correction 
terms to the Newtonian potential coming from an analytic f{R) model, one must consider both the 
second and the fourth order in the perturbation expansion of the metric. Let us consider the perturbed 
metric with respect to a Minkowskian background -f The metric components can be 

developed as follows: 

r) ~ -1 + g^\t, r) + r) 

grr{t,r) ~ l + gil\t,r) 

< ( 11 ) 

geeit.r) = 

= r^sin^6> 

where c = 1, = ct ^ t, and a general spherically-symmetric metric has been considered: 

= g^rdx'^dx'^ = —gQQ{x'^ ,r)dx^‘^ + grr{x^ ■,'r)dr‘^ + r'^dkl ( 12 ) 

where dVt is the solid angle. The next step is to introduce an analytic Taylor expandable f{R) functions 
with respect to a certain value i? = Rq. 

fiR) = E ^^^(R - RoT ^ /o + f,RoR + f,RmR^ + fmR' + ••• (13) 

n\ 

n 

In order to obtain the weak limit, we must introduce Equations (11) and (13) in the field Equations (2) 
and (3) and to expand the equations up to the orders 0(0), 0{2) and 0(4). This approach permits us to 
select the Taylor coefficients in Equation (13). We remember that at zero order 0(0), the field equations 
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give the condition /o = 0, and thus, the solutions at further orders do not depend on this parameter. The 
field equations at (2)-order become: 

/,R0^“ ‘^f,R0gtt!r + ^f,RRoR^? ~ f,R0'’^9tt\r + = 0 


“ ‘^f,R0grl]r + Sf^RRoR^r^ — f,R0^9tt]rr — 0 


( 2 ) 


,(2) 


< ‘^f,R0grl^ — f,R0^~ f,R09tt]r ~ f,R09ri]r + ^f,RRoR^? + ^f,RR0fR 


,(2) 


( 2 ) 


>( 2 ) 


?( 2 ) 

',rr 


= 0 


(14) 


/,r?o^ + Qf,RR0 


‘Igf) + r 


^gfi - rR^^^ + 2g^rilr + rglf. 


2R^? + rR[?r 


( 2 ) 


= 0 


,(2) 


= 0 


As we can see, the fourth equation in the above system (the trace) gives us a differential equation with 
respect to the Ricci scalar, which allows solving the system exactly at 0(2)-order. Then, we obtain: 


( 2 ) _ . <^ 2(^)6 

9tt — Oq --h 


—Lr 


f,Ror 


3L Lr 


+ 


m 


QL'^ 


Lr 


a^^'> - 

' yrr — 


^1 _ ^2(t) Lr + 1 ^_^^ ^ Ssjt) 


f,Ror 


R^^^ = 62{t)- 


^—Lr 


3L Lr 
63(1) e^'' 


6A2 Lr 


(15) 


+ 


2L r 


where L = 


df.RRO 


, f^Ro and f^RRQ. In the limit f{R) —)■ R, we recover the standard Schwarzschild 
solution. Let us stress that L has the dimension of length, the integration constant (5o is dimensionless, 
while 5i{t) has the dimensions of length~^ and (52(f) the dimensions of length~‘^. Since the differential 
equations in the system Equation (14) contain only spatial derivatives, we can fix the functions of time 
5i{t) (i = 1, 2) to arbitrarily-constant values. Furthermore, we can set 5o to zero, because it is not an 
essential additive quantity. If we want a general solution of previous Equation (15), we exclude the 
Yukawa growing mode, obtaining the following relations: 


< 

ds^ = — 

1 

1 


dR + 

1 + 

(52(f) Lr + i ; 


f,Ror 

' 3L Lr 


f,Ror 

3L Lr 


dr'^ _l_ 


(16) 


R = 


52{t)e 


—L r 


where rg = 2MG. Now, we can look for the solution in terms of the gravitational potential from 
Equation (15), obtaining the following relation: 


$ = 


GM ^ (52(f) 


f^Ror ' 6L Lr 

Eet us note that the L parameter is related to the effective mass: 


m = 


= 

L2 j 


2 /, 


RRO 


/. 


RO 


(17) 


(18) 
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and it can be interpreted also as an effective length. Being 1 + — f,m and (5i 

quasi-constant, the Equation (17) becomes: 


6GM 5 


L2 1 + 5 


r \ \jrlvl / —I-\ 

■ (TT+ 


(19) 


Here, the first term is the Newtonian-like part of the gravitational potential to point-like mass, and the 
second term is the Yukawa correction including a scale length, L. If 5 = 0, the Newtonian potential 
is recovered. With this assumption, the new gravitational scale length L could naturally arise and 
accommodate several phenomena ranging from Solar System to cosmological scales. 

3. Constraining/(/?) Gravity Models Using Clusters of Galaxies 

Clusters of galaxies are the largest virialized object in the Universe. They are the intermediate step 
between the galactic and the cosmological scales; thus, any relativistic theory of gravity must be capable 
of correctly describing their physics. Clusters typically contain a number of galaxies ranging from a 
few hundreds to one thousand, grouped in a region of ~2 Mpc, contributing 3% to the total mass of the 
cluster. A more important component is represented by the baryons residing in a hot inter cluster (IC) 
gas. Although IC gas is highly rarefied, the electron number density is Ue ~ 10“^ — 10“^cm“^; it makes 
up 12% of the total mass; and it reaches high temperatures ranging from 10^ to 10® K, becoming a strong 
X-ray source with a luminosity typically ranging Lx ~ 10^® — lO'^^erg/s. 

One of the most promising tools to study clusters of galaxies is the Sunyaev-Zeldovich (SZ) 
effect [106,107]. CMB photons cross clusters of galaxies, and they are scattered off by free electrons 
present in the hot IC gas. This interaction produces secondary temperature fluctuations of the CMB 
power spectrum due to: (I) the thermal motion of the electrons in the gravitational potential well of the 
cluster (it is named Thermal Sunyaev-Zeldovich (TSZ) ; [106]); (2) the kinematic motion of the cluster 
as a whole with respect to the CMB rest frame (it is named Kinetic Sunyaev-Zeldovich KSZ ; [107]). In 
the direction of a cluster (n), the SZ effect is given by: 



( 20 ) 


where dr = aTU^dl is the optical depth, ar is the Thomson cross-section, is the Boltzmann constant, 
rrteC^ is the electron annihilation temperature, c is the speed of light, z/ is the frequency of the observation 
and Vci is the peculiar velocity of the cluster. Tq = 2.725 ± 0.002 K [108] is the current CMB blackbody 
temperature, and g{u) is the frequency dependence of the TSZ effect. In the non-relativistic limit, g[x) = 
a;coth(a;/2) — 4, with x = hu/ksT the reduced frequency. The physical description of the TSZ is 
commonly given by introducing the Comptonization parameter: 



( 21 ) 


where Pe{r) is the electron pressure profile that must be specified to predict the TSZ anisotropies. Using 
X-ray observations and numerical simulations, several cluster profiles (rigTe) have been proposed. 
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Isothermal /3-model fits the X-ray emitting region of the elusters of galaxies well [109,110], with the 
electron density given by: 


?7,e(r) = ne,o 


1 + 


3/J 


( 22 ) 


where is the central electron density and is the core radius. Those parameters, together with 
the electron temperature and the slope (3, are determined from observations. Using the X-ray surface 
brightness of clusters, the slope value ranges in the interval [ 0.6 — 0 . 8 ] [ 111 ]. 

The isothermal /3-model over-predicts the TSZ effect in the outskirts of the cluster of galaxies 
[112]. Recently, a phenomenological parametrization of the electron pressure profile, derived from the 
numerical simulations, has been proposed [113,114]. The functional form of Universal pressure profile 
is: 

Pq 

(c5ooa:)^“[l + 

where x = r/rsoo, and rsoo is the radius at which the mean overdensity of the cluster is 500-times the 
critical density of the Universe at the same redshift. Then, C 500 is the concentration parameter at 
The model parameters were constrained using X-ray data [114,115] or CMB data [116]; their best fit 
values are quoted in Table 1 . 


Table 1. Parameters of universal pressure profile fitted by different groups using 
X-ray [114,115] and CMB data [116]. 


Model 

C 500 


Pa 

la 

Po 

Reference 

Arnaud et al. 2010 

1.177 

1.051 

5.4905 

0.3081 

8.403 h%^ 

[114] 

Sayers et al. 2013 

1.18 

0.86 

3.67 

0.67 

4.29 

[115] 

Planck et al. 2013 

1.81 

1.33 

4.13 

0.31 

6.41 

[116] 


Since TSZ does not depend on redshift in the ACDM model, clusters are a very important laboratory 
to test cosmology (see [117-120] and the references within), even to test the fundamental pillars of the 
Big Bang scenario [121,122]. 

In the last decade, SZ multi-frequency measurements have been reported by several groups: 
the Atacama Cosmology Telescope (ACT) [123-126], the South Pole Telescope (SPT) [127-130] and 
the Planck satellite [1 16,131-133]. Wilkinson Microwave Anisotropy Probe (WMAP) seven years [134] 
and Planck [116,131] data have discussed several apparent discrepancies with ACDM predictions. If 
these discrepancies are due to the physical complexity of clusters of galaxies or are due to the limitation 
of the theoretical modeling [135], it is an open question. 

To study if these discrepancies are due to the theoretical modeling and, at the same time, to 
constrain/rule out ETGs at cluster scales, the pressure profile of a cluster of galaxies has been considered. 
On the one hand, it is interesting to study models that can explain the cluster without requiring any 
dark component [50]. On the other hand, models constructed to mimic the ACDM expansion history 
by replacing the DE with higher order terms in the gravitational Lagrangian must also describe the 
emergence of the ESS and the dynamical properties of clusters [51-53]. 
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3.1. Pressure Profile from Yukawa-Like Gravitational Potential. 


One approach to test ETGs using the TSZ pressure profile is focused on the f{R) Lagrangian that 
is expandable in Taylor’s series (Equation (13) in Section 2.2), where the Yukawa-like correction to the 
Newtonian potential, Equation (19), could be interpreted as the dynamical effect of DM in clusters of 
galaxies [50]. Making the hypothesis that: the gas is in hydrostatic equilibrium within the modified 
gravitational potential well (without any DM term): 


dP{r) 

dr 

the gas follows a polytropic equation of state: 


-p{r) 


d^{r) 

dr 


(24) 


P{r) (X p^{r) 


(25) 


the Equations (19), (24) and (25) can be numerically integrated to compute the pressure profiles by 
closing the system with the equation for the mass conservation: 


dM{r) 

dr 


= 47rp(r) 


(26) 


Thus, the pressure profile will be a function of the two extra gravitational parameters (5, L) and the 
poly tropic index 7. Eet us remark that the density p(r) includes only baryonic matter, without resorting 
to DM particles. 


3.1.1. Data and Results. 


To constrain the model, the Planck (all Planck data are publicly available and can be downloaded 
from http://www.cosmos.esa.int/web/planck) 2013 data and an X-ray clusters catalog [136] have been 
used. In Eigure 1, the predicted Comptonization parameter for the universal pressure profile with 
the parameters given in Table 1 and the = 2/3-model is compared with the /(i?)-model with the 
parameters [5, ^, 7 ] = [—0.98,0.1,1.2]. The model is particularized for the Coma cluster, located 
at redshift z = 0.023, with core radius Tc = 0.25 Mpc, X-ray temperature Tx = 6.48 keV and 
electron density 71^,0 = 3860 m“^. The plot shows that the higher order terms in the gravitational 
/(i?)-Eagrangian could potentially explain the cluster without accounting for DM contribution. 

The SZ emission was measured on the Spectral Matching Independent Component Analysis (SMICA) 
map (SMICA is a foreground-cleaned map; it has been constructed using a component separation method 
by combining the data at all frequencies [137]; the SMICA map has a 5' resolution) at the locations 
of all X-ray selected clusters, and then, it was compared with the model prediction. The temperature 
anisotropies were averaged over rings of width 6 * 500 / 2 , where 6*500 is the angular scale subtended by the 
r 5 oo- The error associated with each data point was computed carrying out 1000 random simulations 
and evaluating the profile for each simulation. The analysis was performed by computing the likelihood 
function log£ = — x ^/2 as: 

X^(p) = Sfj=o(|/(p, Xi) - d{xi))C~^{yfp, Xj) - d{xj)) (27) 

where N is the number of data points, and p = (5,7/, 7). In Equation (27), d{xi) are the data, 
and Cij is the correlation matrix between the average temperature anisotropy on the discs and rings. 
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Two parameterizations have been tested: (A) L = (r^oo, the seale length is different for each; clusters 
scale linearly with rsoo; (B) L has the same value for all of the clusters. The models were constructing 
choosing “physical” fiat priors on the parameters. Looking at the Yukawa-potential in Equation (19), 
if 5 < — 1 , it becomes repulsive, and if <5 = — 1 , it diverges; while the poly tropic index 7 was varied 
in the range corresponding to an isothermal and an adiabatic state of a monoatomic gas, respectively. 
The priors are summarized in Table 2. 



Figure 1. Comptonization parameter for Comacluster (z = 0.023). The pressure profile 
is integrated along the line of sight for: the three universal profiles (dashed, solid and 
dash-dotted lines; the model parameters are quoted in Table 1); /9 = 2/3 model (long dashed 
line); and the f(R) model (red solid line, [5, L, 7] = [—0.98, 0.1,1.2]). 

Table 2. Priors on the parameters when modeling the cluster profile in /(i?)-gravity. 


Parameterization 

S 

7 

L 

C 

(A) 

[-0.99,1.0] 

[1.0,1.6] 

- 

[0.1,4] 

(B) 

[-0.99,1.0] 

[1.0,1.6] 

[0.1,20] 

- 


Figures 2 and 3 show the 2D contours at the 68 % and 95% confidence levels of the marginalized 
likelihoods for both Parameterizations (A) and (B), respectively. Although the contours are opened and 
only upper limits can be obtained, the value 5 = 0 is always excluded at more than the 95% confidence 
level. Therefore, Newtonian potential without DM cannot fit cluster pressure profile, and thus, either DM 
or modified gravity must be the right description of the dynamical properties of the cluster of galaxies. 
The upper limits are summarized in Table 3. 
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Table 3. Upper limits on the model parameters for the eluster profile in /(i?)-gravity. 



68% CL 

95% CL 

68% CL 

95% CL 

6 

<-0.46 

<-0.10 

<-0.43 

<-0.08 

7 

>1.35 

>1.12 

>1.45 

>1.2 

L(orC) 

<2.5 

<3.7 

<12 

<19 



Figure 2. 2D eontours at the 68% (dark green) and 95% (light green) eonfidenee levels of the 
marginalized likelihoods for Parameterization (A). In panel (a), (b), and (e) there are shown 
the 2D eontours for the parameters (C, L), (5, L), and (5, Q, respeetively. Since the contours 
are opened, only upper limits on the parameters can be given. 
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Figure 3. 2D contours from the marginalized likelihoods for Parameterization (B). Contours 
follow the same eonvention of Figure 2. This parameterization also provides opened 
eontours. Therefore, also in this ease, one ean only give upper limits on the parameters 
and ean not distinguish whieh parameterization is the best one. 

3.2. Chameleon Gravity: Hydrostatic and Weak Leasing Mass Profile of Galaxy Cluster. 

A seeond approaeh in using elusters of galaxies to test ETGs is foeused on the ehameleon gravity 
models that introduce a sealar field non-minimally eoupled to the matter eomponents. This eoupling 
gives rise to the fifth foree, whieh is tightly eonstrained at Solar System seales (see Seetion 2.1). 
The main differenee with the other approaeh diseussed above is that the ehameleon gravity only 
eliminates DE, preserving the DM eomponent as the souree of the gravitational field needed to explain 
the emergenee of the ESS. 

Beeause of the sereening meehanism, the effeet of this foree is negligible at small seales, and thus, 
it does not affeet the density distribution of the IC gas in the eores of galaxy elusters, but its effeet may 
not be totally sereened in the outskirts. As a eonsequenee, the gas distribution will be more eompaet 
in the outskirts under the effeet of an additional pressure term that balanees the effeet of the fifth 
force. Thus, the hydrostatie mass of a eluster should be different from the mass obtained from weak 
gravitational lensing that depends only on the distribution of DM along the line of sight and does not 
assume hydrostatie equilibrium. 

The model studied in [51,52] assumes the spherieally-symmetrie distribution of IC gas and DM. 
The IC gas is in hydrostatie equilibrium within the potential well generated by the DM distribution. 
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Departures from hydrostatie equilibrium are parameterized, introdueing a non-thermal term in the 
pressure, while the ehameleon field direetly eontributes to the total amount of mass. The model is 
eompared to leasing observations used to estimate the mass, with X-ray and SZ observations used to 
reeonstruet the gas pressure profile [51-53]. 

Let us start writing down the equation of the hydrostatie equilibrium for the IC gas eomponent: 

1 ^ GM{< r) 

Pgas{r) dr r2 

where M{< r) is the total mass enelosed in a sphere of radius r, pgas is the gas density distribution within 
the sphere and Ptot is the total gas pressure that ean be reeast as: 


Ptotir) = Pth(r) + Pnon-th(?^) + P<i,ir) (29) 

ineluding both thermal (Pth(r)) and non-thermal (Pnon-th(^")) pressure, and the term due to the 
ehameleon field (P<^(r)). Thus, the total mass in Equation (28) is: 


M(< r) = Mth(r) -f M^on-thij) + M^{r) 

By definition, eaeh term ean be expressed as: 


Mth(r) = - 
M'non-th(?^) 

M^(r) ^ - 


dPthir) 
Gpgas(r’) dr 

_ ^ dPnon—thij'^ 

Gpgas{r) dr 

jS d(i){r) 

GMpi dr 


(30) 


(31) 

(32) 

(33) 


M^{r) is the mass assoeiated with the ehameleon field, and Mnon-th represents the departure from the 
hydrostatie equilibrium mass Mth. Introdueing the equation of state of the IC gas. 


Pth{r) = /cngas(r)Tgas(r') 


(34) 


the Equation (31) ean be re-written as 


Mth(r) 


kTgas{r)r / dlnpgas(r) dlnTgas(r’) 

firripG \ dlnr dlnr 


(35) 


where the identity Pgas{r) = prUpUgas^r) has been eonsidered with the mean moleeular weight p and the 
proton mass nip. 

Aeeording to hydrodynamieal simulations of the ACDM model [138,139], the non-thermal pressure 
term ean be modeled as: 


P Jr) - 

4 non—thl' J — / \ 

1 - g{r) 


ngas{r)kTgas{r) 


(36) 


with: 


9 {r) 


ant{l + zf^^ 



M 


200 


3 X IO^^Mq 


riM 


(37) 
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where a^t, j3^t, n^t and nyi are eonstants. Therefore, the non-thermal mass term is: 

M tv, 1 — ^ ^ f ng) 

non-thermal — dr \l - g J ^ ^ 

This term is modeled without considering the chameleon field, and it has been shown to reproduce 
accurately the non-thermal contribution, even in /(i?)-gravity. Nevertheless, a more accurate procedure 
would require modeling this term from hydrodynamical simulations under chameleon gravity. Finally, 
the total mass inferred from hydrostatic equilibrium has to be equal to the the total mass as inferred from 
weak lensing (WL): 

Mtot = ^thermal + ^non-thermal + = MwL (39) 

The DM density distribution is equally well described in the /(i?)-model and in Newtonian gravity 
using the Navarro-Frenk-White (NFW) profile [140]. This profile has the following functional 
form [141]: 


p{r) 


Ps 

r/rg (1 -f r/r^f 


(40) 


and it has been constructed using iV-body simulations of the ACDM model. Thus, the DM mass within 
the radius r is: 


pr 

^^dm(< r) = 47r / drr‘^p{r) 

Jo 


47ip^r^ 


ln(l -f r/vg) 


r/r^ 

1 -f r/rs 


= 3Twl 


(41) 


that is equivalent to the WL mass if the DM component dominates over the baryonic one. The model has 
been tested using only COMA cluster [52] and using a sample of 58 X-ray selected cluster associated 
also with weak lensing measurements [53]. 


3.2.1. Data and Results 


The chameleon model is essentially described by the following parameters (/9, 0oo) in Equations (5) 
and (7). In order to normalize them, those parameter have been replaced by: 


/52 = 


fJ 


1 + / 3 ’ 


^ 00,2 


= 1 — e 


(42) 

(43) 


and constrained by carrying out the Monte Carlo Markov chain (MCMC) method using X-ray, SZ and 
WL observations. Another important parameter for studying clusters of galaxies is the critical radius: 


'^crit — 




(44) 


where ps is the density at this radius. It represents the distance from the DM halo center where the 
screening mechanism acts and it is not possible to distinguish between chameleon gravity and GR [51]. 

In [52], the authors have constrained chameleon gravity only using the Coma cluster (z = 0.0231). 
For this cluster, they implement the MCMC method using the X-ray temperature profile reported by the 
X-ray Multi-Mirror Mission (XMM) -Newton [142] and Suzaku [143] for the inner and outer region. 
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respectively. Then, the X-ray surface brightness profile was taken from XMM-Newton [144] and the SZ 
pressure profile measured by Planck [131]. For the WL counterpart, the measurements were reported in 
[145]. More recently, in [53], the authors have constrained the chameleon model using publicly-available 
weak leasing data provided by the Canada France Hawaii Leasing Survey (CFHTLenS; [146]), and 
X-ray data taken from the XMM-Newton archive. Their sample includes 58 clusters spanning the redshift 
range 2 ; = [0.1 — 1.2], with measured X-ray temperatures = [0.2 — 8] keV. For this sample of clusters, 
they also carried out an MCMC analysis. In Figure 4 are shown the two-dimensional contours for the 
parameters (32 and 0oo,2 as generated by [53]. The red and blue lines are the 95% and 99% confidence 
limits from [52]. Those results represent the best constraints on chameleon gravity using clusters. On the 
one hand, at low values of {3, the departure from GR is too small to be detected with these observational 
errors. On the other hand, GR gravity is recovered within the critical radius rcrit, and large values of (3 
lead to a lower value for (poo that could keep r^rit in the whole cluster region. Although the analysis from 
[53] using the profiles outside of the critical radius of the cluster could better constrain large values of 
(3, the analysis carried out by [52] also used SZ measurements and could discriminate better between 
GR and chameleon gravity at lower values of {3. In both [52,53], it also constrained the term / ro from 
Equation (7). The result provides an upper limit today |/,ro| < 6 x 10“^ at a 95% of confidence level. 


0.8 


0.6 

8 

3 - 

0.4 


0.2 

Figure 4. The confidence contours for the renormalized parameters {(poo, 2 ,132) The 95% and 
99% confidence levels are plotted in light gray and medium gray, respectively. The results 
are from [53]. The 95% and 99% confidence contours from [52] are over-plotted also in red 
and blue, respectively. The vertical line corresponds to |/,ro| < 6 x 10“®. 
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4. N-Body Hydrodynamical Simulations inf(R) Gravity 

All chameleon-like theories of gravity need a speeifie sereening meehanisms (Seetion 2.1). The effeet 
of all sueh meehanisms is refieeted on the non-linear perturbations of the matter distribution. Therefore, 
some of the most optimal tools to provide a full deseription of these meehanisms are modifieations of 
standard N-body algorithms [147-153]. 

The aims were related to the study of the non-linear effeets on the emergenee of the LSS, sueh 
as the matter power speetrum [148,150,154], the redshift-spaee distortions [55] and the physieal and 
statistieal properties of the eold DM halos and voids [56,140,155-160]. However, the same seales 
and struetures that eould provide a powerful tool for testing the theory of gravitation are also seriously 
affeeted by non-linear astrophysieal proeesses [161-166]. In the few last years, a very powerful eode 
to investigate the eombined effeets of astrophysieal proeesses and alternative theories of gravity has 
been developed [57]. The eode is named Modified Gravity (MG)-GADGET. It was used to study the 
Hu-Sawieki model [103], pointing out the existenee of an important degeneraey between modified 
gravity models and the uneertainties in the baryonie physies partieularly, related to AGN feedbaek 
[57]. Afterwards, the eode was eombined with a massive neutrinos algorithm [167], to study the eosmie 
degeneraey between modified gravity and massive neutrinos. It was found that the deviations from the 
matter power speetrum in the ACDM model are, at most, of the order of 10% at z = 0 taking f{R) 
models with = — 1 x 10“"^, and the total neutrinos mass = 0.4 eV [168]. Using again the 

Hu-Sawieki parameterization, the effeets of the model on the IC gas and its physieal properties have 
been studied. It was found that the DM veloeity dispersions in the f{R) model are ~4/3 larger than in 
ACDM, when low mass halos are eonsidered. Nevertheless, the mass of the objeets that are affeeted by 
the ehameleon field depend on the value of the baekground field: for \f^Ro\ = 10“^ the field is never 
sereened, nor for massive elusters; for |/,ijo| = 10“®, the ehameleon field is sereened for mass above 
~1O^^'^M0; and for |/,ijo| = 10“®, the Hu-Sawieki model does not produee any effeet in the range 
of mass explored in [169]. Finally, MD-GADGET and the Hu-Sawieki model have been also used to 
study the synthetie Eyman-a absorption speetra using simulations that inelude star formation and eooling 
effeets. The diserepaneies with the simulations based on the ACDM model were at most 10% with the 
baekground field fixed at \f^Ro\ = 10“'^ [58]. 

All tests earned out using N-body simulations have provided several indieations about the presenee of 
a degeneraey between the effeets of the eomplex baryonie physies (highly non-linear) and the effeets of 
the modifieations of the theory of gravitation. Even with the high aeeuraey of future data, without a full 
understanding of the non-linear astrophysieal proeesses, it would be not possible to break the degeneraey 
and to eonstrain or rule out ETGs. Therefore, there is the need to fully explore this degeneraey with 
N-body simulations to know the level of bias introdueed. 

5. Constraining the Expansion History of the Universe in f(R) Gravity 

Cosmologieal models ean be tested looking at baekground evolution and also at the growth of 
the eosmie struetures. Testing the f{R) model at the baekground level allows finding the so-ealled 
“realistie” models that have the radiation, matter and DE eras very elose to the ACDM ones. Although 
many realistie f{R) models have been eonstrueted, it should be eonsidered that different models with a 
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comparable late-time aeeelerated expansion and similar baekground densities’ evolution eould differ in 
the growing of the matter density perturbations [170]. For example, it is well known that the effects at 
the baekground level of eonsidering a eollisional matter component filling the Universe are negligible. 
Nevertheless, non-eollisional and eollisional matter affeet in different ways the evolution of matter 
perturbations, opening the possibility of distinguishing them [171]. Another problem is represented 
by the need to break the degeneraeies between the evolving DE models and modified gravity in order to 
confirm/rule out models [172]. Thus, the growth faetor (GF) data become extremely important to test 
modified gravity and to find their signatures. However, the available data have less accuraey than the 
other datasets, sueh as SNela, Baryonie Aeoustic Oseillation (BAG) and thus, it is not possible to 
fully eonstrain /(/?)-models. This scenario will totally ehange with the fortheoming Euclid observations. 
That dataset will have an unpreeedented aeeuraey, and it will allow eonstraining deviations from the 
standard model at a few pereent level [54]. 

Let us start by the flat Friedman-Lemaitre-Robertson-Walker metrie (FLRW), 

= —dt^ + a^{t)dx^ (45) 


where a{t) is the rate expansion of the Universe. The Ricci scalar is related to the Hubble funetion, 

H = by the following expression: 

a{t) 


R = 6 (2H^ + id) 


(46) 


and using the non-relativistic matter and radiation approximation, the following eonservation laws 
are satisfied: 


Pm + Pm — 0 (47) 

-f 3Hp:y = 0 (48) 

Using Equations (2), (3), (45) and (46), the modified Friedman equations are obtained [173]: 

= n\p^ + p^) + i(/,^ii - f\R)) - 3Hf,n (49) 

—‘^f,RH = K^{pm + -p-y) + f,R — Hf^Ji (50) 


To compute the the evolution of the matter density perturbations, dm = one ean mainly follow 
two different approaeh: 

• Sub-horizon approximation: The evolution of the matter density perturbations can be obtained 
solving the following equation [174]: 

dm + 2i7dm - 47rGefT(a, k)p^5^ = 0 (51) 


where k is the comoving wave number, and Gefr(a, k) is the effeetive gravitational eonstant. The explieit 
forms of Geff depend on the ETGs model. In the ease of f{R) gravity, it is given by [175]: 


G (fc^/g^) {f^Ru/f^n) 

f^Rl 1 + 3 (/c7«^) if,RR/f,R) . 


Geff(a, k) 


(52) 
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By definition, the growth factor fg is: 

^ din5^ 

^ din a 

thus, Equation (51) becomes: 

dfg . [ ^ , r\ r _ 3 STrGeffPm 

2 3H^ 


(53) 


(54) 


Since Equation (51) does not give raise to analytical solutions, it is customary to introduce an efficient 
parameterization of the matter perturbations. The most common one is the so-called 7 -parameterization: 


fo 




(55) 


jACDM obtained for the value 7 ~ 0.545. Therefore, any departures from this value could indicate the 
need to use a different model. More in general, the 7 function could evolve with the redshift in the case 
of modified gravity models [176,178,179]: 

7 (^)= 70 + 717 — (56) 

This approach has been widely used to constrain the f{R) gravity, pointing out that: the current data 
have not enough constraining power to discriminate between different modified gravity models and the 
ACDM model; the so-called growth index is the most efficient way to provide constraints on the model 
parameters; and the 7 -parameterization represents the best way to measure deviation from the standard 
model [176,177]. 

• Dynamical variables: Starting from Equations (49) and (50), it is possible to reach a 
model-independent description of the evolution of the matter density perturbations [54,180,181]. 
Defining the following dimensionless variables: 


Xi = 

X2 = 


f,R 

Hf,R 

fjR) 


X 3 = 

X4 = 


R 


the density parameters and the effective equation of state of the system can be written as: 


(57) 

(58) 

(59) 

(60) 


III 

(61) 

^ ^ Pm ^ 

- 0 rr9 x ^2 X 3 X 4 

j,R 

(62) 

Utx = xi+ X 2 + X 3 

(63) 

'^eff = -^(23^3 - 1) 

(64) 
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The background evolution is given by: 


x\ = —1 — Xs — 3 X 2 + x\— X 1 X 3 + X 4 

X 1 X 3 


x '2 = 


X 3 = 


m 

X1X3 

m 


- X2{2X3 - 4 - Xi) 


2 x 3 (X 3 - 2 ) 


X 4 = — 2 X 3 X 4 — X 1 X 4 


(65) 

( 66 ) 

(67) 

( 68 ) 


where the prime indicates the derivative with respect to the scale factor d/d\na, and m is a parameter 
related to the theory that indicates how much a specific model is close to the ACDM (or how much it 
deviates from it). This parameter is defined as: 


d f,R _ Rf,RR 

dlnR 


(69) 


and by definition, one also finds the following identity: 


d\nf _ _R£r _ ^ 
dlnR f X 2 


(70) 


Equations (65) to (70) represent a closed system that can be integrated numerically, obtaining the 
evolution of the background densities perturbations. We just need to choose an /(i?)-model to compute 
the two parameters m and r that we need to integrate the system. From the analysis of the critical points, 
we have several conditions that an f{R) model has to respect to be cosmologically viable [54,180]. 
In general, only models with m > 0 and very close to the ACDM model (f{R) = R — 2A) are 
cosmologically viable. Thus, the equations for the matter density perturbations are [181,182]: 


C + (^ 2:3 - ^Xi^ 5 ;, - ^(1 - Xi - X2 - X3)5„ 


— - 6 + 3xi - 3 X 4 - 3xi(x3 - 1) > df^R 


X 


+ 3(-2xi + X 3 - l)5/,ij + 3(5/, 


R 


df,R + (1 - 2xi + X3)5f^R 

^ 2 x 3 

— — 2x3 H- 

l_X 5 m 


xi(x 3 + 1 ) - X 4 + x\ 


df, 


R 


= (1 - Xi - X 2 - X3)6m - Xi(5^ 


(71) 


(72) 


where (5 /,_r = df^n/f^R and X 5 = aH satisfies: 


X5 = (X3 - 1)X5. 


(73) 


The evolution at all scales of the growth of the matter density perturbation, dm, could be obtained 
via numerical integration. However, a more straightforward approach is to assume that the growth 
rate, fg, is a function of time, but not of scale [183-187], and to use the 7 -parameterization in 
Equations (55) and (56). 
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The detection of possible deviation from ACDM model has been studied fitting the 7 -parameterization 
to the forthcoming Euclid dataset [54]. Assuming the ACDM model as the reference model, it has been 
found that: the forthcoming data will be able to constrain the parameter 7 and w with an accuracy of 4% 
and 2% if they do not depend on redshift (standard evolution) and to clearly distinguish the ACDM from 
alternative scenarios at more than 2(j; see Figure 5. 


1.0 


0.5 


n 0.0 

-0.5 


- 1.0 


Figure 5. The la and 2a marginalized contours for the parameters 70 and 71 in the 
7 -parameterization. Shown is the reference case (shaded yellow regions), with the optimistic 
error bars (green long-dashed ellipses) and the pessimistic ones (black dotted ellipses). 

Red circle represent the ACDM model (7 = 0.545), while triangles represent the f{R) 
model [54]. 

6. Testing Gravity Using the Cosmic Microwave Background Data 

The CMB temperature anisotropies play a key rule to constrain the cosmological model. They 
are mostly generated at the last-scattering surface, and thus, they provide a way to probe the early 
Universe [13,17,133]. However, the lack of a full-consistent theoretical framework and the fine tuning 
problem related to the cosmological constant. A, have encouraged studying wide classes of alternative 
DE or modified gravity models. The most difficult challenge of modern cosmology is to use current and 
forthcoming datasets to discriminate between those alternative visions of the Universe [54,188-191]. 
Many of those models mimic the ACDM at background evolution, but differ in the perturbations, and 
they could produce several effects on the CMB power spectrum. As an example, different models 
could have different expansion histories that shift the location of the peaks [192], and could affect 
the gravitational potential at late time, changing the integrated Sachs-Wolfe (ISW) effect [193,194]. 
Moreover, modification of the gravitational theory should be reflected in modifications of the lensing 
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potential [195,196], the growth of the struetures [197-200] and also the amplitude of the primordial 
B-modes [201,202]. 

Modified gravity models affeet the evolution of the Universe at both the baekground and the 
perturbation [203-206]. There are mainly two different approaehes used to eonstrain alternative models 
of gravity: the parametrized modified gravity approaeh and the effeetive field theory. 

The parametrized modified gravity approaeh eonstruets funetions probing the geometry of spaee-time 
and the growth of perturbations. Starting from a spatially-flat FLRW metrie, the perturbed line element 
in the Newtonian gauge is given by: 

ds^ = a(r)2[-(l + 2^)dT^ + (1 - 2^)dx^] (74) 


where r is the eonformal time. 4/ and $ are the two sealar gravitational potentials that are a funetion of 
the seale ^ = ^(k, a) and $ = <h(/c, a) and are related to the energy-momentum tensor, through 
Einstein’s equations. Those potentials ean be modeled to fix all degrees of freedom in order to mateh 
observational data. Moreover, the differenee of those two potentials <h — 4/, whieh is the anisotropie 
stress, has to be equal to zero in GR. Any departure from this value would indieate a departure from 
GR [207,208]. One of the most eommon parameterization is implemented in the Modifieation of 
Growth with Code for Anisotropies in the Mierowave Baekground (the eode is publiely available at 
http://www.sfu.ea/ aha25/MGCAMB.html) (MGCAMB) [209,210], and it eonsiders the following: 


= —4:7iGa‘^^{k,a)pA (75) 

^ = ^[k,a) (76) 

to parameterize the Poisson and the anisotropie stress equations using two seale and time-dependent 
funetions n{k, a) and 7 (fc, a). Those two funetions are determined by ehoosing the partieular model. 
For example, some elasses of f{R) models lead to [211]: 


/i(/c, a) 

7(fc,a) 


1 -f /SiA^ 

1 -f Af k'^a^ 

1 + (^ 2^2 

1 -f A| 


(77) 

(78) 


where the parameters fdi are dimensionless and represent the eouplings, and the parameters A* are lengths. 
Those funetional forms have been used in early studies to foreeast errors of the parameters for the 
Dark Energy Survey (www.darkenergysurvey.org/) (DBS, [212]) and Large Synoptie Survey Teleseope 
(http://www.lsst.org/lsst/) (LSST, [213]) surveys [210]. The main limitation of this approaeh is the use 
of the quasi-statie regime assuming that the seales at whieh one is interested are still linear and smaller 
than the horizon. Furthermore, the time derivatives are negleeted. 

Effeetive field theory (EFT) [214,215] deseribes the whole range of the sealar field theories. The 
Lagrangian preserves the isotropy and homogeneity of the Universe at the baekground, providing a more 
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general approach without assuming the quasi-static regime, but adding much more free parameters that 
must be constrained. The /(i?)-models are a sub-class of those theories. The action of EFT reads: 


^ = 


■ a c{T)og 


[1 + ^(r)] R + A(r) 

+ MM - M!(T)2a%«'SKj; 

- + ^UdMsg«>sR‘-^^ 





Sm [x*) 9fJ.u\ 


(79) 


where is its spatial perturbation, Kj^ is the extrinsic curvature, uiq is the bare (reduced) 

Planck mass and Sm is the matter part of the action, including all fluid components, such as 
baryons, cold DM, radiation, and neutrinos, but it does not include DE. The nine time-dependent 
functions [216] {Q, c, A, Mf, M|, M|, M|, m^}, have to be fixed to specify the theory. The EFT has 

been implemented in the publicly-available EFTCAMBcode (http://www.lorentz.leidenuniv.nl/hu/codes/. 
Version 1.1, October 2014.) [59,217], that numerically solves both the background and perturbation 
equations. The code has been used to constrain massive neutrinos in f{R) gravity [59]. 

Planck collaboration has used in both approaches to fit f{R) models to investigate their effect 
combining early time datasets, such as CMB, and cosmological datasets at much lower redshift [14]. 
In general, f{R) models require setting the initial condition: 


lim 

R^oo 


m 

R 


0 


(80) 


and the boundary condition: 


B{z) 


fRR HR 
I + IrH-H^ 


( 81 ) 


High values of Bq, its present value, change the ISW effect and the CMB lensing 
potential [211,218,219], Marchini2013. Their analysis has shown the presence of a degeneracy between 
the optical depth r and Bq that can be broken adding other probes on the structure formation, such as 
weak lensing, CMB lensing or redshift-space distortions (see Figure 6). 
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Figure 6. The 68% and 95% contour plots for the two parameters, {Logio(Bo), r}. There 
is a degeneracy between the two parameters for Planck temperature power spectrum (TT) 

+ BSH (the combination of BAO, SNIa, and Hq datasets). Adding lensing will break the 
degeneracy between the two. Here, Planck indicates Planck TT. 

The same study carried out using MGCAMB code led to the same results within the uncertainties, 
improving the previous constraint on Bq [220]. 

Nevertheless, both EFTCAMB and MGCAMB assume a specific background cosmology. MGCAMB 
code assumes for the background evolution the ACDM or teCDM, while EFTCAMB assumes a fixed 
background evolution, and the CMB power spectrum is calculated for an f{R) model re-constructed by 
the fixed expansion history. A more recent code, named f{R) Code for Anisotropies in the Microwave 
Background (FRCAMB) [221], allows reconstructing the CMB power spectrum for any f{R) model 
by specifying the first and second derivatives of the Eagrangian. Although the results obtained are 
consistent with the one from EFTCAMB and MGCAMB, the code has the advantage of being designed 
for any f{R) model at both the background and perturbation level. 

7. Discussion and Future Perspectives 

In this review, we have outlined the possibilities offered by the current and the forthcoming 
astrophysical and cosmological dataset to constrain or rule out the /(i?)-models. It is mandatory to 
search for the answer to one of the fundamental questions in modem cosmology: is GR the effective 
theory of gravitation? This question comes from the evidence that GR is not enough to fully explain 
the cosmological evolution of the Universe and the emergence of the clustered structures, and it needs 
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the addition of two unknown components. The difficulties in identifying the DM and DE components 
at a fundamental level have opened the possibility that the answer resides in modifying the theory of 
gravity. Nevertheless, having an alternative demands testing it study the physical phenomena at much 
smaller scales than the cosmological ones. To this effect, much analysis has been carried out to ensure 
that modified gravity could recover the tight constraints at the Solar System scale. Models that survive 
these probes have been used to test galactic and extragalactic scales. Clusters of galaxies provide an 
excellent laboratory to probe the different features of the different modified gravity models. The variety 
of the phenomenology related to the cluster physics allows one to use SZ, X-ray and WL observations 
in order to constrain /(i?)-models. Many efforts have been also devoted to quantitatively describing 
the physical effects of alternative cosmological scenarios. Those efforts have required developing new 
algorithms for hydrodynamical N-body simulations, to study the degeneracy between baryonic physics 
and modifications of gravity and also new Boltzmann codes to predict the CMB power spectrum and to 
take advantage of the Planck data to constrain modified gravity models. 

Nevertheless, a self-consistent picture has not been reached yet. The current datasets do not have 
enough statistical power to discriminate between modified gravity models and standard cosmology, 
and when the needed precision is present, several degeneracies have been found. Next generation 
experiments, such as the LSST [213] and the Wide-Field Infrared Survey Telescope (WFIRST) 
(http://wfirst.gsfc.nasa.gov/) [222], will undertake large imaging surveys of the sky, while other surveys, 
such as Euclid [223] and the Big-Baryon Oscillation Spectroscopic Survey (BigBOSS) [224], will make 
spectroscopic measurements of galaxies. The combination of all forthcoming datasets will hopefully 
put strong constraints on the evolution of the Universe and the emergence of the ESS, allowing us to 
discriminate between the concordance ACDM model and its alternatives. 
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